A GPS Navigator for the Free Energy Landscape. 
Used to Find the Chirality-Switching Salt 
Concentration of DNA 

Joshua T. Berryman* and Tanja Schilling 

University of Luxembourg, Luxembourg 



E-mail: josh.berryman@uni.lu 



Abstract 

Sets of free energy differences are useful for finding the equilibria of chemical reactions, 
while absolute free energies have little physical meaning, however finding the relative free 
energy by subtraction of absolute free energies is a valuable strategy in certain important cases. 
We present calculations of absolute free energies of biomolecules, using a combination of the 
well-known Einstein Molecule method (for treating the solute) with a conceptually related 
liquid method of recent genesis (for treating the solvent and counterions). The approach is 
based on thermodynamic integration from a detailed atomistic model to one which is simplified 
but analytically solvable, thereby giving the absolute free energy as that of the tractable model 
plus a correction term found numerically. An example calculation giving the free energy with 
respect to salt concentration for the B- and Z-isomers of duplex DNA in explicit solvent and 
counterions is presented. 
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Introduction 

Overview: Recent improvements of atomistic forcefields for nucleic acids and advances in a 
particular family of thermodynamic integration techniques'^ together encourage a return to what 
is becoming a canonical model system to benchmark computational methods for biomolecules and 
poly electrolytes: the chirality- switching B to Z isomerisation of DNA. 

A new addition to the progression of free energy methods beginning from the "Einstein Crystal" 
of Frenkel & Ladd,^ followed by the "Einstein Molecule",^ the "Confinement Method"^ and a 
recently introduced analytical reference-liquid modeP^ is presented here and demonstrated for 
B and Z DNA. In the absence of an existing umbrella term for this philosophy of free energy 
calculation, it is referred to here as Thermodynamic Integration to Analytical Solution (TIAS) 
because it relies on performing a Thermodynamic Integration (TI) calculation with a realistic (but 
intractable) model at one endpoint of the integration path and an analytically solvable (AS) model 
at the other. The principal technical advance within the field of TIAS to be presented here is in the 
adaptation of a recently developed reference liquid model 5 8 such that it can be used for triangular 
water in a molecular dynamics simulation, and successful combination of this adapted liquid model 
with the existing Einstein Molecule (EM) model for the solute. 

For the 'realistic', but intractable, end of the integration path, the AMBER forcefield and sim- 
ulation code were used. The details of this quantitative all-atom treatment (with Coulomb, bonded 
and van der Waals interactions) are given in the methods. 

The results of the example calculation are found to be good, although scope for further devel- 
opment is acknowledged. 

B-Z DNA Isomerisation: The B-Z isomerisation of duplex DNA is a gross structural transi- 
tion, from a right-handed double helix (the canonical B-form) to a left-handed double helix (the 
Z-form, which occurs in-vitro at high salt concentration^^ or in-vivo with the assistance of DNA- 
binding proteins 11 and/or with negative supercoiling 12 ). Occurrence of Z-DNA in mammalian 
cells has been implicated as a causative factor for certain cancers. 13 In low salt the Z-form is dis- 
favoured relative to the B-form mainly for reasons of electrostatics and solvation; it becomes free- 
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energetically favourable in total when a high salt concentration screens the repulsion between the 
charged backbone phosphate groups, also making partial solvent-exposure of the base pairs more 
favourable (Fig.[T]). In formal terms a first order phase transition is not possible in this system, as 
it is one dimensional with finite interaction ranges. However the unwinding which is required to 
change the handedness of the double helix presents a formidable energetic barrier which must be 
overcome for the transition to take place.LL^LS 




Figure 1: B and Z isomers of [d(CG)g]2 DNA, after 12ns equilibration in 1 Molar NaCl. The 
backbone phosphate groups are closer together in the more compact Z isomer, and the base-stacks 
have greater solvent exposure. 

The salt concentration for coexistence of these two forms is known from circular dichroism 
measurements to be roughly 2.5M NaCl (depending on sequence, temperature, chain length and 
methy lation) J 9 1 1 ° l 1 ^ The availability of this experimental data, particularly for short dodecamer 
sequences, has attracted several simulation studies at atomistic levels of detail. Definitive ther- 
modynamic characterisation with respect to salt concentration has remained just over the horizon 
for some years, however recent improvements have been made: 3DRISM has been qualitatively 
successful in finding the coexistence point (estimated at 1M); and targeted molecular dynamics 
simulations have had at least some success in describing the transition path * 14 ! 15 ! although this is 
difficult to compare directly against experiment. 

Thermodynamic Integration to Analytical Solution: This established and powerful family 
of thermodynamic integration (TI) techniques has recently been extended to permit application to 
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liquids J^EIS] In brief, the method consists of defining for a given substance a reference model which 
has an analytically tractable expression for its free energy. The Helmholtz free energy Aq of the 
system under the realistic (but intractable) Hamiltonian is then expressed: 



A =Ai- 




where A\ and M{ refer to the free energy and Hamiltonian of the reference system; and /o() 
and are mixing functions used to control the speed with which ^ is introduced and is 
removed on the interval A = [0, 1]. 

The scheme for solids or macromolecules, due to Frenkel and Ladd 6 with modifications by 
Vega, 7 is known as the "Einstein molecule method" (EM) or as "the confinement method". The 
basic idea is to introduce a reference model in which the particles do not interact with each other, 
such that the reference partition function factorizes and the reference free energy can easily be 
computed exactly. In the simplest such model, the particles are simply coupled to lattice sites by 
harmonic springs. 

If one is interested in locating a phase transition, one can, of course, use eq. 1 (or a similar 
expression) in order to compute the free energy difference between two phases directly and avoid 
the "detour" via absolute free energies. However, eq. 1 only holds if the integration path does not 
cross a first order phase transition. The capability gain of absolute free energy calculation (as op- 
posed to performing TI over single or multi-legged free energy cycles only between non-tractable 
Hamiltonians) is thus that states which are separated by a first order phase transition can then be 
compared. This presents a much simpler alternative to the practice of carrying out TI calculations 
which integrate between phases by following a complex system- specific path designed to give con- 
sistently pseudo-critical behaviour.^^ Additionally, with one endpoint of the integration fixed for 
all future calculations, some of the "black art" of path design which has previously given TI on 
biomolecules a reputation for difficulty and expense should become susceptible to being automated 
away as integration paths for calculations on radically different systems come more to resemble 
one another. 
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Einstein Molecule & Cage-Swapping Liquid Model: The EM model used for the DNA 
molecules is a simple noninteracting-particle treatment in which each atom rests in a harmonic 
well, and has been thoroughly discussed elsewhere.^ The reference model used here for the water 
and counterions is still quite new, however, and has been slightly modified from the most recent 
appearance in the literature,^ for use in conjunction with MD rather than MC dynamics. 

The EM model of harmonic wells cannot be used as a reference for the liquid state, because 
in the liquid each particle is ultimately free to move throughout the volume; hence each particle 
would have to sample a well of infinite width when approaching X = 1. Instead we use an attractive 
potential of short range (Fig. [2], here called 'CS', with reference to its 'cage-swapping' dynamical 
behaviour which is explained below). In addition to motion within the potential by normal MD, 
an MC move was incorporated allowing pairs of particles to exchange the wells to which they 
are attracted. This 'Cage-Swapping' move reflects a third, factorial, term in the general partition 
function of liquids, related to the indistinguishability of particles. In effect each time a pair of 
particles exchanges wells they are exchanging identities, thus respecting the statistics of the liquid 
state. At the same time, the cage-swapping move serves to enhance the sampling by allowing 
particles to exchange wells without needing to cross the saddle point which separates them; for a 
detailed explanation see Schilling & Schmidt 

Modifications made for the MD/MC CS model used here relative to previous pure-MC treat- 
ments were to use MD timesteps instead of the 'displacement' and 'relocation' MC moves; and to 
add a harmonic region near to the well-bottoms; such that the wells were harmonic out to 1 A, then 
constant-force out to 5A, then flat. The small harmonic region was necessary for stability of the 
MD. 

A natural objection to the use of a swap move in conjunction with MD simulation is to ask if it is 
required that the MC sampling should converge before each MD timestep takes place; mercifully it 
has been shown in other work that any combination of MD and MC steps which would individually 
provide Boltzmann sampling of the degrees of freedom to which they are applied will in total 
provide Boltzmann sampling of the combined degrees of freedom.^ 
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Figure 2: A simple harmonic well was used for the EM potential, while a constant-force well with 
a cutoff at long range and a harmonic region at very short range was used for the CS potential of 
the fluid part of the system. 

Motivation of TI Method: In a TI calculation, the free energy landscape of the system grad- 
ually morphs from one shape to another, in a way that is not yet easy to characterise. From an 
intuitive basis, we can state qualitatively that for effective TI four sometimes-conflicting criteria 
are important: 

(a) To minimise the error from the numerical integration, the integration path should be as 
smooth as possible; meaning that the higher-order derivatives of each generalised force 
4jf{X)M > should be as small as possible. 

(b) The shapes of landscapes to be mixed should be similar to each other; otherwise the gen- 
eralised force due to Mq, for instance, could sporadically take extremly large values when 
/o(A) is small and J£f is primarily controlling the dynamics. 

(c) To minimise the convergence time per integration point, exploration of the mixed landscapes 
should be as fast as possible. 

(d) Because the variance of the generalised force scales with the gradient of the mixing function 
/(A), to efficiently distribute sampling between integration points, the mixing functions 
should be as close to linear as is consistent with the other criteria. 

To meet criterion (a), it is first required to avoid discontinuities in the generalised forces by 
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making sure that the integration does not cross any first order phase transitions. The basis of TIAS 
is a successful treatment of this issue by defining tractable models which match a given phase in 
terms of the symmetries expressed by their partition functions. Definition of the model partition 
functions is discussed in the supplementary information. 

Criterion (b) is more difficult. Given that the landscape under J£f is in fact topologically 
different to that under Mq, because particles can pass through one another under M\, the question 
of appropriately changing between landscapes is brought sharply into focus. In order to treat cases 
where a cavity in the free energy landscape (e.g. a region of infinite potential energy, such as due 
to the Lennard- Jones short range repulsion) must be introduced, starting from a flat potential, it 
is valuable to gradually increase the maximum energy associated with that region to some high 
but finite value. The region is only then blocked completely, by introducing the infinite potential 
inside the region which is already effectively forbidden by the large finite potential. This strategy, 
which was attempted here by use of the 'Guide Hamiltonian' discussed below, should ensure that 
regions of phase space which are populated for a given A are also populated at neighbouring values 
of A, while still allowing topological change of the landscape to occur. In the absence of such a 
technique, the generalised force can diverge, with an infinitesimal change in A giving an infinite 
change in < f(X)J^ >. 

A further, partial, fix for problems with respect to criterion (b) is provided by choice of the 
mixing function used to introduce or remove a given Hamiltonian. Careful work has been done 
on this in the context of alchemical TI involving Lennard- Jones and Coulomb forces,^ which was 
broadly followed here: the essence of the strategy is to use polynomial mixing functions such that 
as the ability of a given Hamiltonian Jtf? to direct the system into its own minima drops with /(A), 
the generalised force Jr/(A) also drops. The choice of mixing function must be balanced 
between criteria (b) and (d). 

Criterion (c) is addressed here by defining J£f, the CS Hamiltonian, with few and low energetic 
barriers, such that exploration becomes almost freely diffusive over the course of the integration. 
It has been suggested that it might also sometimes be valuable to retain barriers in the landscape 
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as these can serve to reduce the effective dimensionality of the configurational space and actually 
accelerate exploration, 21 however this effect was not observed here. 

Criterion (d) is addressed here by choosing a relatively low-order mixing function compared to 
some of those which have been trialled in other work.^1 

Guide Hamiltonian: To manage the changes of the topology of the free energy landscape 
over the integration, a third Hamiltonian, J^ gU ide^ was introduced in addition to Mq and M\, with 
a mixing function such that it would be zero at both endpoints of the path. The weak short-range 
interparticle repulsion provided by Jtf'guide controlled the topology in the sense that it was used to 
first introduce gently sloping energy barriers around those states which would later take on infinite 
energy becoming topological 'cavities' in the free energy landscape. 

This use of a guide Hamiltonian, introduced here, is an alternative to the practice of adding a X- 
dependence to the functional form of the Lennard- Jones (LJ) potential such that it is finite over all 
configurations for A > 0, known as 'softcoring'.^Softcoring of the LJ typically follows mixing- 
out of the Coulomb interaction (which also has a discontinuity at zero separation), in a multi-step 
procedure which increases computational expense. This process is also undesirable because it has 
the potential to introduce first order phase transitions: the phase diagram of water, for instance, is 
dramatically altered when electrostatics are removed. 

The use of an M'guide has the advantages of solving at a stroke discontinuities of both LJ and 
Coulomb forces, of being extremely simple to implement and to parameterise, and of introducing 
no extra code into the complex and highly optimised non-bonded force and Ewald sum routines of 
the molecular dynamics program. The short-range repulsion is handled using a neighbour list and 
therefore adds only a small (and linearly scaling) cost to the calculation. 

The guide Hamiltonian was defined very simply as a quadratic repulsive potential between all 
heavy atoms, where the cutoff was defined as the radius giving the minimum of the van der 
Waals interaction for the given atom pair ij, multiplied by 0.88. 
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^uide = £ £ / MIN[^ ; -r ;j ,0] \ 2 
Mixing Function Design 

Numerical smoothness of the integration over A depends critically on the choice of mixing func- 
tions used to turn the Hamiltonians Mq, 3tf{ and ^ K ^e °ff an d on - The mixing functions were 
defined as polynomials: 



/o(A) = (1-A) 4 (3) 
= A 2 (4) 
fpadefr) = ^A 2 (l-A) 4 (5) 

The choice of polynomial mixing functions (which has been advised with order 4 or greater for 
alchemical transformations involving Lennard- Jones atoms 20 ) gives generalised forces j^fo(X)J% 
and jifi which are zero for A = 1 and A = respectively, so that (for example) two parti- 

cles can pass through each other without causing a singularity in the generalised force with respect 
to Mq- 

Formally, the use of a polynomial mixing function is enough to stabilise the integration, mixing 
the two topologically distinct landscapes without introduction of singularities. Unfortunately the 
multiplication of very large by very small numbers is numerically abhorrent, therefore the guide 
Hamiltonian was also required. 

In order to allow ^f gu id e to take effect before removing the non- smooth Hamiltonian 3%q, the 
mixing function /o(), was rescaled such that the endpoint of the mixing would be at A = X g instead 
of at A = 1 (with A g chosen arbitrarily as 0.5). The actual function used in place of /o() was 
therefore: go(X) = MIN[/o(^),0]. The three mixing functions are shown in Fig.|3j 
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Figure 3 : The function go (A ) controlled the AMBER Hamiltonian Mq, f\ (A ) controlled the analyt- 
ically tractable reference Hamiltonian J£f and f gu ide(k) controlled the 'path' Hamiltonian J^ gu ide- 



Results 

Free Energies: Free energies were separately calculated for Z and B conformations over a range 
of salt concentrations (Fig. [4]). The crossover between the B and Z regimes occurred somewhere 
near the 2.8M NaCl found in the experiments of Pohl and Jovin.^ As far as the authors are aware, 
this is the most accurate numerical calculation of the coexistence position to date. 
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Figure 4: Free energy differences between the B and Z conformations with respect to salt concen- 
tration. Coexistence points found by other groups using 3D-RISM (MYH)^ and experimentally 
(PJ)E are indicated with arrows. 



An interesting feature of the results which would probably not be apparent in experiments 
is that away from the coexistence point, the relative stability of the B-form appears to decrease 
slightly for very low salt. Although this behaviour lies within the errorbars it could be a hint 
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of the Z-B-Z pattern of stability with respect to salt concentration which has been observed for 
methylated 

It was found that, while Z-DNA is strongly metastable at low salt concentrations, the more 
labile B-DNA structure tended to melt quickly at 2.75M NaCl and above. For this reason data 
points at these very high salt concentrations were discarded. 

Integration Path: The character of the integration path can be seen by examining the gener- 
alised force due to each of the three Hamiltonians with respect to A (Fig. [5]). The smoothness of 
these curves (apart from at A =0.5, where the integration-out of the AMBER Hamiltonian is com- 
pleted) indicates that no first-order phase transitions were crossed during the integration path, and 
also serves to justify the use of a basic Simpson's rule integration scheme; carried out separately 
over the three terms of the generalised force and over the [0,0.5] and the [0.5,1.0] intervals. 

1000 
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Figure 5: Generalised forces due to the AMBER Hamiltonian (J%), the EM+CS Hamiltonian 
(M\) and the 'guide' Hamiltonian used to prevent particle overlaps during the mid-section of the 
integration (J4? g ). The forces due to the two main Hamiltonians were so large that the difference 
between B and Z is not visible on this plot. The guide Hamiltonian (right axis) was much smaller 
than the main Hamiltonians (left & right axes). 
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Methods 

AMBER Setup: Molar NaCl concentration was defined as the ratio of Na + ions to H2O molecules 
multiplied by the molarity of water, 55.55. was defined using the AMBER99 forcefield with the 
Barcelona corrections to nucleic acid parameters,^ the Joung-Cheatham ion parameters^ and the 
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TIP3P water model.^ Water molecules were kept rigid using the SHAKE algorithm. Temperature 
was held at 300K using a Langevin thermostat with a coupling of 0. lps -1 . The molecular dynamics 
timestep was 2fs. 

Start configurations of DNA in the B and Z conformations were prepared using the NAB molec- 
ular building tool.^Counterions were added near to the DNA using Xleap so as to neutralise back- 
bone charge, and further counterions were then added at random to bring the salt concentration up 
to the required molarity assuming 9441 water molecules. Water molecules were then added from a 
pre-equilibrated water box, deleting those which overlapped with an ion or DNA atom, until 9441 
were present. The resulting configuration was energy-minimised and then equilibrated at 300K 
and 1 atmosphere in the NPT ensemble for 2ns. The average volume of each system was taken 
over the interval 1 to 2 ns, and the system box sizes in each configuration were set to these average 
volumes, so that the TI calculations could be carried out in the NVT rather than NPT ensemble. 

Analytical Model Parameter Choices: In the reference model of the system, solute atoms 
were placed in harmonic wells of spring constant 5kg TA . Counterions and water oxygen atoms 
were assigned to CS wells, which had a spring constant of 1.112kgTA~ 2 up to r\ = 1A, then a 
constant force of 1.1 1kg TA -1 out to r2 = 5 A. The 'cage- swapping' reassignment of CS wells was 
carried out using a Metropolis MC algorithm. For each liquid molecule present in the system, 5 
MC well-swapping attempts were made every timestep. 

Molecular images were prepared using VMD.^ 

Performance and Convergence 

Calculations were carried out using four 2.26 GHz processor cores per integration point (therefore 
68 per measurement). Calculations required on average 26.3 hours per lOOps block giving a total 
of 210,000 core-hours for the entire calculation. 

To estimate the convergence of the individual generalised force measurements at each integra- 
tion point, a spectral analysis of the time series of generalised forces was carried out using the 
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CODA library of R functions, in order to identify the decorrelation times and effective numbers 
of independent samples present. Estimated Standard Errors of the mean (ESE) are given as the root 
of the variance divided by the effective number of independent samples. The ESE at each value of 
A (Fig. [6^b)) provides a description of the efficiency of the sampling in each of the 'mixed' free 
energy landscapes which were visited. 




0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 

( a ) A W A 

Figure 6: (a) Variance and Decorrelation Time of Generalised Force. The rate of convergence 
of the estimate of the generalised force at a given A was determined by variance (left axis) and by 
the decorrelation time T (shown on the right axis). Numbered labels indicate peaks in the variance 
or in T which are discussed individually in the text, (b) ESE of the Generalised Force. The 

ESE per 500ps was much higher for A < 0.5, therefore longer runs of 1600ps were made for this 
section. 



The radically different sampling efficiencies at the different integration points derive from sep- 
arately addressable causes, the signatures of which are individually labelled in Fig. [6} and which 
are discussed in relation to the criteria for good TI design which were presented earlier: 

(a) Smoothness of Generalised Force w.r.t. A : The smoothness criterion was met quite well 
(Fig.©. 

(b) Small Variance of Generalised Force: When the potential corresponding to a given term 
of the generalised force is very different to the potential which is controlling the dynamics, 
then fluctuations in that term of the force can become very large, the peak labelled '4' on 
Fig. [6] shows this effect. 
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Peak '4' is lower than in first attempts, thanks to the use of a polynomial mixing function, and 
also greatly diminished by the use of the guide Hamiltonian. Peak '4' might also be further 
reduced by the use of different Hamiltonians which more closely mirror the dynamics of the 
real system. 

(c) Fast Exploration: The peaks labelled '2' and '3' show slow exploration due to roughness of 
the potential energy landscape: sampling of the Boltzmann distribution for atomistic mod- 
els of biomolecules is famously difficult, due to the many traps and barriers which exist. 
Sampling also becomes slower in the CS Hamiltonian when it is fully turned on, due to the 
depth of the wells becoming greater than the thermal energy. Peak '5' shows a slowing of 
exploration due to decoupling of particles under the CS Hamiltonian: when the particles no 
longer interact, decorrelation over time must rely purely on the Langevin thermostat rather 
than being assisted by the chaotic properties of the many-body Hamiltonians. 

To address point (c) there is no easy escape from the need to sample the atomistic model at 
one end of the integration ('2'). At the other end of the integration ('3'), the trapping due 
to the CS potential can be decreased by using shallower wells, but at the cost of making the 
configurational space which is explored more different to the space of the atomistic system. 
Although the aspect indicated by peak '5' seems to be of relatively minor importance, it 
should be possible to address it by altering the coupling parameter of the Langevin thermostat 
for different values of A, or by using an M{ which has explicit interparticle forces while 
remaining tractable. 

(d) Near-linear Mixing: Peak '1' arises because the gradient of the mixing function g{X) is 
large at small A, and this large value scales not only the generalised force but its variance 
- this represents a tradeoff accepted in the process of mixing function selection, it stands in 
balance against peak '4'. 

To address this need to trade between the criteria (b) and (d) it seems that the best hope is 
to improve the guide or endpoint Hamiltonians used such that the need for a polynomial mixing 
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function is removed. 

Discussion 

This calculation represents the first quantitative numerical estimate of the salt concentration at B-Z 
DNA coexistence, at least as far as the authors are aware. As recently as 2010 it was noted that 
such a calculation should be prohibitively difficult "because the problem requires the sampling of 
an extraordinarily large configuration space, including water and ions, to obtain the free energy 
difference";^ that the result is quantitatively correct demonstrates solid improvements in the tech- 
niques of free energy estimation suitable for biomolecular complexes, over and above the normal 
progress due to improved hardware. It is anticipated that the software and parameter sets devel- 
oped can be reused widely and can be effective for any combination of biomolecule, counterions 
and solvent (and also that there is still substantial room for further development of them). It is 
intended to release the software and parameter sets as a downloadable library on Dr Berryman's 
website, for easy linking against the AMBER simulation program or other packages. 

The accuracy of the final result serves as a validation for use of the combination of the Joung- 
Cheatham ion parameters and the parm99Bsc forcefield (both of which are relatively new) with 
TIP3P water for simulation of DNA in the fairly unusual conditions of very high salt concentration. 

The work here has made minor alterations to the 'Cage-Swapping' model for absolute free en- 
ergy calculation of fluids and combined it with a standard Einstein Molecule method for the solute. 
The principal novelty was in the scale and complexity of the system treated, but the introduction 
of a 'guide' Hamiltonian to control the path of the thermodynamic integration while not altering 
behaviour at the endpoints is also novel as far as the authors are aware. 

There is too much active research in thermodynamic integration and free-energy perturbation 
methods to give a comprehensive list of ideas which could influence further work. A brief sug- 
gestion is that absolute free energy methods could be well suited to mapping of phase diagrams 
of complex fluids. There are a large number of biomolecular and soft-matter systems which ex- 
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hibit rich phase behaviour that have not yet been explored. On the methodological side, recent 
advances in nonequilibrium TlJ^in correlated sampling^ and in Hamiltonian-exchange 31 all of- 
fer increases in computational efficiency for future TIAS calculations. It has been suggested that 
an efficiency gain can be made by doing away with the intermediate stages of the TI calculation 
entirely, instead carrying out a type of importance sampling over and such that those areas 
of configurational space which overlap between them receive enhanced attention.^ 

The most obvious step in further development, and one which should be complementary to 
most other avenues of research, is extension and refinement of the analytical models which serve 
as endpoints for the calculation. Ideally, the mixed energy landscapes generated by an analytically 
tractable model should be as similar as possible to 'softened' versions of the initial Hamiltonian. 
This possibility should serve as a call to arms for any theorists who may believe that they can 
swiftly assemble a tractable model for some system of interest. 

A supporting document is provided, giving derivations for the absolute free energies of the 



CS and EM Hamiltonians. This material is available free of charge via the Internet at http : 
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